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Abstract 



> 

^N . We study the use of the Evans Nonequihbrium Molecular Dynamics (NEMD) 

heat flow algorithm for the computation of the heat conductivity in one- 

On . dimensional lattices. For the well-known Fermi-Pasta-Ulam (FPU) model, it 

J>^. is shown that when the heat field strength is greater than a certain critical 

Q . value (which depends on the system size) solitons can be generated in molec- 

,^ . ular dynamics simulations starting from random initial conditions. Such soli- 

tons are stable and travel with supersonic speeds. For smaller heat fields, no 
solitons are generated in the molecular dynamics simulations; the heat con- 
ductivity obtained via the NEMD algorithm increases monotonically with the 
size of the system. 
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I. INTRODUCTION 



Heat conduction in one- dimensional (ID) lattices has attracted much research interest. 
Surprisingly, it has been found that many ID lattices do not obey the Fourier's law |[l|]- | T7| : 



the thermal conductivity is divergent in the thermodynamic limit. For harmonically coupled 
oscillators, it was rigorously shown that the thermal conductivity A is proportional to the 
number of oscillators A^ |I[]. Such a divergence is founded in the existence of extended waves 
(phonons) freely traveling (and carrying thermal energy) along the lattice without attenua- 
tion. In later studies [0,0, impurities or defects in the chain were taken into account, since it 
was anticipated that phonon waves could be damped by the scattering processes due to such 
defects, thus possibly removing the A^ divergence of A. However, it was demonstrated for iso- 
topically disordered harmonic chains that the heat conductivity still diverged at a somewhat 
slower rate (A ~ N^/'^) @,0]. Another way of trying to achieve normal heat conduction in 
one-dimensional lattices is by invoking anharmonicity [Q: here nonlinearity makes it possible 
for phonons to interact among themselves thus impeding their free propagation. However, 



Lepri et al [0 have found that even strong nonlinearity and chaotic behavior is insufficient 
to ensure the existence of normal heat conduction. In the well-known Fermi-Pasta-Ulam 
(FPU) /5-model they found a power-law divergence of the thermal conductivity A oc A^''' for 
7 K. 0.4. This power-law divergence was qualitatively attributed to the long-time tail of the 
heat flux autocorrelation function, whose time integral gives the thermal conductivity of the 



system []13 



Previous studies of heat conductivity have used a straightforward simulation method 



0- |[T^: In molecular dynamics (MD) simulations two heat reservoirs with high and low 
temperatures Th and Tl respectively are located on each side of the lattice. The average 
heat flux and the internal temperature gradient are measured, with the thermal conductivity 
being the ratio of these two quantities. However, there are a number of disadvantages with 
this approach. In particular, the system is spatially inhomogeneous and one cannot define an 
intrinsic temperature T for the system due to the large temperature gradient. Consequently 



it is impossible to obtain the T dependence of the heat conductivity. In addition, problems 
associated with the use of the Nose-Hoover thermostat for boundary particles have been 



discussed in Ref. |TH] 



Recently Maeda and Munakata [^ proposed a homogeneous nonequilibrium molecular- 
dynamics (NEMD) method based on the Evans heat flow algorithm [p!^-pD|. However,the 



system size used in Ref. |T0| was too small (32 particles) to allow one to study the behavior 
of the NEMD algorithm in the thermodynamic limit. 

In this paper we present a detailed study of the Evans NEMD heat flow algorithm for 
one-dimensional (ID) lattices. We demonstrate that when the heat field is sufficiently large, 
well-defined solitary waves (solitons) can be generated in simulations with random initial 
conditions. These soliton objects travel at a supersonic speeds, and they also appear in 
the corresponding Hamiltonian systems. When a soliton is generated, the normal process 
of homogeneous heat conduction is destroyed, as heat is transferred through the chain via 
the highly localized solitary wave. This results in a sharp increase in the effective thermal 
conductivity of the system. Due to this instability, progressively smaller fields are required, 
as the system size increases, to observe the linear regime of the thermal conductivity and 
thereby carry out the extrapolation of the thermal conductivity to zero field strength. 

The paper is organized as follows. In section 2 we describe the NEMD equations of 
motion for one-dimensional lattices with nearest neighbor interactions, and we point out the 
possibility of solitary wave solutions in the system. In section 3 we carry out the nonequi- 
librium heat fiow simulations on the well-known Fermi-Pasta-Ulam (FPU) model. We show 
that for the heat field strengths greater than a certain critical value, which depends on 
the system size and temperature, solitons can be generated from random initial conditions. 
For smaller heat fields, no solitons can be generated in molecular dynamics simulations with 
random initial conditions; in this case, the heat conductivity can be obtained via the NEMD 
algorithm. Some concluding remarks are presented in Section 4. 



II. NEMD EQUATIONS AND SOLITARY WAVE SOLUTIONS 

We consider a ID system of N particles located along the x axis with lattice constant 
a = 1. Each particle is allowed to move in the y direction, perpendicular to the x-axis, and 
we denote by g, the displacement of the i-th particle and by pi the corresponding momentum. 
The Hamiltonian of the system is expressed as 



iV 



j=i 



1 



^ = E ^Pl + U{q^+l-q^) , (1) 



where rrii is the mass of the i-th particle, and U represents the interparticle interaction 
potential. For example, the FPU /?— model is exemplified by U{r) = ^r'^ + ^r^, r = g^+i — gj. 
The Evans thermal conductivity algorithm for general N-particle systems of fluid has 
been described in detail |]T8|,|T9|. The basis of this algorithm is the Green-Kubo relation 
for the thermal conductivity. For a ID system defined by the Hamiltonian model (|I]), the 
thermal conductivity coefficient is given by 

A = }^J^^ £ dt{Ut)UO))eg, (2) 

where ks is the Boltzmann's constant, T is the absolute temperature of the system, and 
L = Na is the system size (The lattice spatial constant a is taken to be 1 in this paper). 
The notation (■ ■ ■)eg denotes an equilibrium ensemble average. The heat flux vector, J^, is 
given by 

Ut) = -^ T.^iU'{Q^^l - *) + U'{q. - g._i)]. (3) 

In the Evans NEMD algorithm the N particle system is coupled to a "heat field" F^. The 
coupling is defined in such a way that the energy dissipation is proportional to Jx^e, i.e., 
dH/dt=LJxFe, and that the adiabatic incompressibility of phase space (AIT) condition is 



satisfied |T^ . The thermal conductivity coefficient can then be found from the ratio of the 



heat fiux to the applied heat field: 



A = lim hm i^^^. (4) 



Fe^O t^oo TFe 
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Here {Jx{t)) is, in principle, a nonequilibrium ensemble average; but in practice, it can be 
replaced by a time average of Jx{t) if the nonequilibrium steady state is unique. 
The equations of motion which satisfy the above conditions are 

qi = Pi/rrii, (5) 

Pi = Fi + FeDi - api, (6) 

where Fi = U'{qi+i — qi) — U'{qi — qi-i) is the total force on particle i due to the nearest 

neighbor interaction, 

1 1 ^ 

A = -l^[U'{q^+l - q^) + f/'(g. - g.-i)] + ^ E[^'(^.+i - ^.) + U'{q, - q,^i)l (7) 

and a is the thermostat multiplier. For a Gaussian thermostat, a is 

1 ^ »• 1 ^ 

« = ^^E^-(^^ + ^eA), K, = -Y.p]/m- (8) 

This Gaussian thermostat ensures that the system's kinetic energy is fixed at the constant 
value Kq. We note that if Fe 7^ 0, the total force exerted on the system by the heat field is 
exactly zero, i.e., FeX^A = from Eq.(|^). Thus the total momentum, J2Pi{t), will remain 
zero for t > provided that its initial value vanishes. 

We have simulated the above described NEMD equations (il)-(i) for a variety of inter- 
particle interaction potentials, and we observe well-defined solitary wave excitations. Before 
describing in detail our numerical results we first analyze why such solitary waves can occur. 

For simplicity, we set rrii = 1 for all i. From Eqs.(§) and @ we obtain 

qi = Fi + FeDi - aqi. (9) 

Introducing a new variable Qi = qi — qi-i, leads to 

Qi = Fi- F,„i + Fe(A - A-i) - aQi 

= f/'(Q.+i) - 2U'{Q,) + f/'(Q._i) - iFe[f/'(Q,+i) - f/'(g.-i)] - aQ, (10) 

This equation, together with Eq.(|^) for the definition of a, forms a set of lattice equations 
for the variable Qi. Note that if Fe and a are set equal to zero, then Eq. (p!oD becomes a 



lattice system which can support stable supersonic solitary waves for many types of nearest 
neighbor interaction potentials U (see, e.g., Refs. [pT|-p6|). But exact analytical solutions 
for the solitary waves are not available except for some rather special potential functions U. 
For example, in the Toda lattice, U{r) = (6/a)[exp(— ar) + ar — 1], the system is integrable 
and exact analytical soliton solutions are known pT| . 

To analyze Eq.(|lO|) we can use the so-called quasi-continuum approximation techniques 
of Refs. [^ll -|26|| to reduce Eq. ([To|) to a partial differential equation: 



d^ d_ 



Q 






dx"^ 12 dx^ 



d 1 93 
+ 



dx 6 dx'^ ' 



U\Q), xg[0,L] 



(11) 



which can be considered as a generalized Boussinesq equation [^ . Here due to the periodic 
boundary conditions used for the system (|^J|), i.e., Qn+i = Qi, The solution of Eq.([TT|) must 
satisfy 



/ Q(x)dx = 0. 
Jo 



(12) 



We assume that Q{x, t) = Qs{x — Vt) + Qo, where Qq is a constant background and Qs is a 
localized soliton profile which vanishes ds, z = x — Vt -^ ±00. Then from Eq.(|T2|) we have 



Qo = -T [ Qs{x-Vt)dx. 
L Jo 



(13) 



Substituting Q{x, t) = Qs(x~Vt) +Qo into Eq . ([TT|) we find that the solitary wave profile 
Qs{z) = Qs{x — Vt) should satisfy 



V 



aV — 

dz^ dz 



Qs 



92 1 d^ 
+ 



dz"^ 12 dz^ 



fjI^I'' 



' dz 6 dz^ ' 



U'{Qs + Qo). 



(14) 



But it does not seem possible to find an analytical solution for this partial differential 
equation for the general case of Qo 7^ 0, Fg 7^ and a 7^ 0. 

However, we note that if the system is large (A^ ^ 1) then Qq must be very small 
according to Eq. (p^) . Moreover, assuming that F^ and a are small parameters, we can set 
Qq = 0, Fe = 0, and a = in Eq.(|T^). After such a drastic simplification, it is possible 
to find approximate solitary wave solutions analytically for various types of interparticle 
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potentials U{Q) pll-|2^. Here we consider a particular example, the FPU potential U{Q) 



^Q"^ + \PQ^- According to Refs. pTl-p^, the solitary wave profile Q{z) can be approximated 

by 



V 13 cosh[2 J3(y2 _ 1)^] 



Once such an approximate soliton solution is obtained then we can go back to the original 
lattice variables {qi,Pi)- Notice that Qi = qi — qi-i, thus 

g* = go + E Qi ~ 90 + / [Qo + Qs{x - Vt)]dx. (16) 

Consequently, the contribution to the kinetic energy from the soliton is 



Ksoi = 2 Ep? = 2 E ^.— 2 loo "^^ [L ^^^^"^ ~ ^^^H = 2^' loo ^'^'^'^' 



(17) 



73/3 

In numerical simulation we observe that when a soliton is generated, small amplitude phonon 
waves give a negligible contribution to the system's kinetic energy. In such a case, we have 



73/? 2 ' ^ ^ 

where T is the absolute temperature of the system which is fixed by the Gaussian thermostat 
(Here the Boltzmann constant ks is set to be 1). For any values of A^ and T, Eq.(0) has 
a unique solution V > \^ which is the (supersonic) speed of the solitary wave. Moreover, 
according to Eq.(|TB|), the soliton speed increases with the system's total kinetic energy. 
In the next section, we present numerical simulation results to confirm these analytical 
predictions. 

III. SIMULATION RESULTS 



Following Ref. 0] we apply the NEMD heat fiow algorithm [Eqs.(|,|)] to the FPU /3- 



model for which Uppuif^) = t^/2 + /?r^/4, and the parameter f3 is taken to be 1 without loss 



of generality. Periodic boundary conditions are always used. Unless indicated otherwise, the 
initial conditions for qi and Pi are obtained by a random number generator. The equations 
of motion are integrated using a fourth-order operator-splitting integrator which conserves 



the system's kinetic energy [|2^. The time stepsize is 5t = 0.005 and the total simulation 
time is between 10^ to 10^ units for each trajectory. 

The first feature of note is that for a given temperature T and particle number A^, 
stable solitary waves can be generated during simulations (which start from random initial 
conditions) if the heat field strength F^ is greater than a certain critical value Fcr- The 
solitary wave travels in the direction of heat flow (to the right in Fig.l) with a supersonic 
speed {Vs > 1). When the soliton is generated the normal process of homogeneous heat 
conduction is destroyed and the heat flux increases drastically (Fig. 2). In such a case, heat 
is transported in the form of a highly localized energy pulse carried by the soliton, and the 
average value of heat flux is nearly independent of Fg. 

We flnd that the soliton's velocity increases with temperature and system size (Fig.3), 
but it is nearly independent of the heat fleld strength Ff,{< 1.0). For example, for A^ = 
100, T = 10, the soliton's velocities are found to be about 4.7, 4.7 and 4.8 for Ff, = 0.01, 0.1, 
and 1.0, respectively. This is in good qualitative agreement with that predicted by equation 



(|T8|). However, it is found that the analytical result ([TsD overestimates the soliton's velocity. 
For instance, for A^ = 200, T = 1, Eq. ([T8|) gives V = 5.586; in the numerical simulation we 
found that the soliton's velocity is about 3.2. The discrepancy is mainly due to the fact that 
Eq.(p!8D is obtained by using a crude approximation of the soliton solution for the lattice 
system (0). Similar to the case of Hamiltonian lattices pT|-p6|, here we also flnd that the 
soliton's amplitude increases with its velocity. 

Although the spontaneous generation of solitary waves (from random initial conditions) 
is observed in the NEMD simulations only when the heat fleld is strong enough, such waves, 
once generated, continue to exist in the system after the heat fleld and thermostat are 
switched off (Fig. 4). The reason for such a behavior is that the soliton is an inherent 
excitation in the FPU /3-model, as explained in the previous section. 
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When the heat field strength Fe is smaller than the critical value no soliton can be 
generated from random initial conditions (See Fig.5). In this case, the time-averaged heat 
fiux {Jx{t)) can be measured (See Fig.6), and the conductivity, limt^oo {Jx(t))/(TFe), can be 
calculated. In Fig. 7 the heat conductivity obtained through the NEMD simulations is plotted 
(Error bars are estimated to be within 10% at most). The NEMD heat conductivity increases 
with the system size. When Fe ^ the conductivity converges to a finite value, which, 
according to the generalized Green-Kubo realtion |[l^j2^, should equal the conductivity 
obtained through Eq.(|^). We have tested this convergence for the FPU system of A^ = 100 
particles. In Fig. 8 we plot the function 

m = -^^ I dT{UT)UO))e,, (19) 

where the ensemble average is obtained by using ten independent trajectories of the length 
10^ units in time. It is clear that when t — > oo, A approaches to a value around 93, which 
is in good agreement with the heat conductivity obtained through NEMD algorithm (See 
Fig.7). 

Finally, we investigate how the critical field strength Fcr for generating solitons depends 
on the particle number N and the system temperature T. We found that Fcr increases with 
T. For example, at A^ = 100 the critical value Fcr is around 0.0054 and 0.0085 for T = 1 
and T = 10, respectively. On the other hand, we found that Fcr decreases monotonically 
with A^ as shown in Fig. 9. For a system of 10000 particles, the critical field is as small as 
0.0005 (accuracy is within ±10~^). This means that an extremely small field has to be used 
in order to observe the linear regime (with no solitons) of heat conduction. However, when 
the heat field is too small the noise-signal ratio [when using Eq.(^)] would become too large 
and thus the efficiency of the NEMD algorithm will be drastically reduced. 

IV. CONCLUDING REMARKS 

In conclusion, we have shown that the Evans NEMD heat fiow algorithm, which was 
designed originally for computing thermal conductivity in liquids, can generate solitons 



when it is applied to ID lattices. In the well-known FPU model, we have shown that when 
the heat field strength is greater than a certain critical value a soliton can be generated 
from random initial conditions. Such a soliton is stable and it travels with a supersonic 
speed which is determined by the system size and temperature. Because of this instability, 
progressively smaller fields have to be used (as the system size increases) to observe the linear 
regime of the thermal conductivity and thereby carry out the extrapolation to zero field. 
This greatly reduces the efficiency with which the algorithm can be used to compute the 
thermal conductivity of large ID lattices. Nevertheless, for small systems, we have found 
that the heat conductivity increases with the size of the system, which is in qualitative 
agreement with the previous finding |1l2|-p!^. 

The present study has been primarily focused on FPU /5-model, but we have checked 
numerically that similar phenomena also exist for other types of ID lattices with distinct 
interparticle interaction potentials (e.g., the Toda potential and Morse potential) and even 
for diatomic lattices. In particular, the spontaneous generation of solitons can be observed 
not only for the nonequilibrium heat fiow systems with Gaussian thermostat but also with 
the Nose-Hoover thermostat and an isoenergetic thermostat. The mechanism underlying 
the observed chaos-soliton transition above the critical field strength still remains to be 



identified, particularly in terms of Lyapunov spectra-shift and phase space contraction fl^ 
induced by the nonequilibrium heat fiow algorithms. 
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FIGURES 
FIG. 1. The evolution of Qi = qij^i — qi, showing the generation of a sohtary wave in a NEMD 

simulation of heat flow in the FPU model. The heat field strength is F^ = 0.006, the system size is 

N = 100 and the simulation temperature is T=l. Note that due to periodic boundary conditions 

the soliton enters into the left end of the lattice whenever it leaves the right end. 

FIG. 2. The instantaneous heat flux in NEMD simulations of heat flow in the FPU model, with 
Fe = 0.006 (solid line) and Fg = 0.01 (dashed), showing a drastic increase of heat flux when a 
soliton is generated about time t = 1000 in each case. This is in contrast to the situation when no 
soliton is generated for Fe = 0.002 (dotted line). Model parameters are A^ = 100, T = 1. 

FIG. 3. The velocity of the solitons in the FPU model with 100 particles (circles) and 200 
particles (squares), respectively. The heat field strength is F^ = 0.01. Lines are for guidance only. 

FIG. 4. Propagation of a soliton in the Hamiltonian FPU lattice without heat field (-Fg = 0.0) 
and thermostat (a = 0). The initial conditions are taken from the last output of Fig.l. N = 100. 

FIG. 5. The same as in Fig.l but for F^ = 0.002, showing that no solitons are generated. 

FIG. 6. Time-averaged heat flux in NEMD simulations of heat flow in the FPU model with 
Fe = 0.002, N = 100, T = 1. 

FIG. 7. The heat conductivity obtained from the NEMD simulation of the FPU lattice with 
N = 100 particles (circles), N = 300 (squares), and N = 400 (triangles), respectively. Simulation 
temperature is T=l. Lines are for guidance only. 

FIG. 8. The time-dependent heat conductivity obtained from the Green-Kubo relation. 
N = 100, and T = 1. 

FIG. 9. The critical field strength for generating solitons, as a function of the particle number. 
The simulation temperature is T = 1. The line is for guidance only. 
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Fig.1, Zhang, Isbister, and Evans 
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Fig. 2, Zhang, Isbister, and Evans 
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Fig.3 Zhang, Isbister, and Evans 
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Fig. 4 Zhang, Isbister, and Evans 
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Fig5. Zhang, Isbister, Evans 
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Fig 6, Zhang, Isbister, and Evans 
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Fig. 7, Zhang, Isbister, Evans 
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